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ABSTRACT 

We present first results on polarization swings in optical emission of blazars obtained 
by RoboPol, a monitoring programme of an unbiased sample of gamma-ray bright 
blazars specially designed for effective detection of such events. A possible connection 
of polarization swing events with periods of high activity in gamma rays is investigated 
using the data set obtained during the first season of operation. It was found that 
the brightest gamma-ray flares tend to be located closer in time to rotation events, 
which may be an indication of two separate mechanisms responsible for the rotations. 
Blazars with detected rotations during non-rotating periods have significantly larger 
amplitude and faster variations of polarization angle than blazars without rotations. 
Our simulations show that the full set of observed rotations is not a likely outcome 
(probability < 1.5 x 10“^) of a random walk of the polarization vector simulated by 
a multicell model. Furthermore, it is highly unlikely 5 x 10“®) that none of our 
rotations is physically connected with an increase in gamma-ray activity. 
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1 INTRODUCTION 

Blazars are active galactic nuclei whose jets are oriented 
close to our line of sight, so that we observe high relativis¬ 
tic beaming of their non-thermal emission and large ampli¬ 
tude variability at all wavelengths. The low-frequency emis¬ 
sion is dominated by synchrotron radiation, and hence is 
highly polarized. The exact polarization fraction and di¬ 
rection depend on the structure of the magnetic field in 
the emitting region, and on the number of emitting re- 
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gions along the line of sight. The polarization direction 
(in the simple case of a single dominant emission region) 
traces (and is perpendicular to) the direction of the pro¬ 
jected magnetic field on the plane of the sky. Already from 
early optical observations, it has been known that polariza¬ 
tion parameters of blazars are variable on daily time-scales 
(Kinman et al. 1966). In general, both flux density and po¬ 
larization exhibit an erratic variability (Angel & Stockman 
1980; Uemura et al. 2010; Ikejiri et al. 2011), which could 
be interpreted as a random walk (Moore et al. 1982). How¬ 
ever, in some cases the electric vector position angle (EVPA) 
of the polarized emission displays long, smooth and mono- 
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Table 1. Selection criteria for the gamma-ray—loud and the control sample. 


Property 

Gamma-ray-loud sample 

Gontrol sample 

2FGL 

included 

not included 

2FGL F(E > 100 MeV) 

> 10“® cm“^ 

- 

2FGL source class 

agu, bzb, or bzq 

- 

Galactic latitude |6| 

> 10° 

- 

Elevation (Elv) constraints^ 

Elvmax > 40° for at least 90 consecutive 
nights in the window June — November 

Elvmax ^ 40° for at least 90 consecutive 
nights in the window April - November 

R magnitude^ 

< 17.5 

< 17.5 

CGRaBS/15GHz OVRO monitoring 

no constraints 

included 

OVRO 15 GHz mean flux density 

no constraints 

> 0.060 Jy 

OVRO 15 GHz intrinsic modulation in¬ 
dex, m 

no constraints 

> 0.05 


^Refers to elevation during Skinakas dark hours 

^Average value between archival value and measured during preliminary RoboPol observations (when applicable) 


tonic rotations which have been observed in the optical since 
the 1980s (Kikuchi et al. 1988). A number of mechanisms 
have been proposed for the interpretation of such events, in¬ 
cluding: stochastic variations of turbulent magnetic fields, 
a shock travelling through a non-axisymmetric magnetic 
field (Konigl & Choudhuri 1985), polarized flares in the ac¬ 
cretion disc (Sillanpaa et al. 1993), two-component models 
consisting of two independent sources of polarized emission 
(Bjornsson 1982), and jet bending (Abdo et al. 2010a). 

Blazars represent the most common class of known 
gamma-ray sources (Nolan et al. 2012; Acero et al. 2015). 
Despite the recent progress in the field, many questions con¬ 
cerning the high-energy emission produced by blazars are 
still under debate. For instance, it is unclear where the 
gamma-ray emitting site is located: within the broad-line 
region (e.g. Blandford & Levinson 1995; Poutanen & Stern 
2010) or well downstream in the jet (e.g. Marscher et al. 
2008; Agudo et al. 2011). 

Recent work showed that at least some large EVPA 
swings can be associated with gamma-ray flares (e.g. 
Abdo et al. 2010a; Larionov et al. 2013b) and therefore can 
possibly provide some insight on the physics of high-energy 
activity. Although such events have triggered an increasing 
interest in polarimetric monitoring of gamma-ray blazars, 
efforts in this direction have been based on selected cases 
comprising statistically biased samples. As a result, a sig¬ 
nificant amount of invaluable polarimetric data sets for a 
large number of sources has been gathered. However, this set 
cannot be used for statistically rigorous population studies 
and, in particular, the investigation of a possible correlation 
between gamma-ray flares and optical EVPA rotations. The 
RoboPol programme (King et al. 2014; Pavlidou et al. 2014) 
has been designed to provide a data set of rotation events in 
an unbiased sample of blazars, appropriate for such studies. 

In this paper, we analyse EVPA rotations detected by 
RoboPol during the first observing season between 2013 July 
and November. After a brief description of observing and 
reduction techniques in Sec. 2, we estimate the frequency of 


EVPA rotations in blazars and list their properties in Sec. 3. 
A Monte Carlo simulation is performed in Sec. 4 in order 
to determine whether the EVPA rotations can be produced 
by random walk processes. In Sec. 5 we study the possible 
connection between an increased activity in the gamma-ray 
band and EVPA swings. Our findings are summarized in 
Sec. 6. 

2 OBSERVATIONS AND DATA REDUCTION 
2.1 Our sample 

A unique feature of the RoboPol programme is that it is 
monitoring a sample which has been selected on the basis of 
strict, bias-free and objective criteria (for detailed discussion 
on the sample construction, see Pavlidou et al. 2014). The 
sample consists of three distinct groups. 

(i) The main (“gamma-ray-loud”) sample is an unbi¬ 
ased subset of a statistically complete flux-limited sample 
of blazars from the second Eermi-LAT source catalogue 
(Nolan et al. 2012). Specifically, we selected all the sources 
in the 2FGL catalogue classified as BL Lacertae objects 
(bzb). Flat Spectrum Radio Quasars (bzq), or active galaxy 
of uncertain type (agu). Applying the selection criteria listed 
in Table 1, we constructed a gamma-ray flux-limited “parent 
sample”. Application of the visibility constraints and field- 
quality cuts resulted in an unbiased subsample of 83 sources, 
among which we randomly selected 62 sources. 

(ii) A “control” sample of 15 “gamma-ray-quiet” sources. 
It constitutes an unbiased subset of a statistically complete 
sample of blazars. It has been drawn from the CGRaBS 
catalogue (Healey et al. 2008) applying the selection criteria 
listed in Table 1. 

(iii) 24 additional sources chosen on the basis of their 
variability characteristics or their presence either in the F- 
GAMMA programme sample or in TeV catalogues. 
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Figure 1. Season length and median cadence for the first season 
data. Broken lines limit areas where rotations at rates of 10, 15 
and 20 deg d“^ can be detected (see Sec. 3.3 for details). Only 
objects with At-median < 20 were left for the more detailed view. 

Although here we present the polarization swings detected 
in all monitored sources during the first RoboPol observing 
season, the statistical analysis in this paper is based only on 
sources from the group (i). 


2.2 Optical observations 

All photometric and polarimetric measurements were done 
at the 1.3-m telescope of Skinakas observatory^ using 
RoboPol, a polarimeter specifically built for the project 
(King et al. 2014). The RoboPol instrument contains a fixed 
set of two Wollaston prisms and half-wave plates, which 
splits each incident ray into four rays with polarization 
plane rotated 45° with respect to each other. Measuring 
relative intensities in pairs of the rays for each object in 
the 13 arcmin x 13 arcmin field, we obtain the fractional 
Stokes parameters q = {R ~ I 2 )/{R + R) = Q/l and 
u = {R — Ii)/{R + Ia) = U/I. Stokes parameter I is cal¬ 
culated as a sum of intensities of all four spots. Since the 
polarization parameters are measured simultaneously, we 
avoid unmeasurable errors caused by the sky changes be¬ 
tween measurements and imperfect alignment of rotating 
optical elements. 

The data presented in this paper were taken with the 
i?-band filter. Magnitudes were calculated using calibrated 
field stars either found in the literature or presented in PTF 
(Palomar Transient Factory) i?-band catalogue (Ofek et al. 
2012) or USNO-Bl.O catalogue (Monet et al. 2003), depend¬ 
ing on availability. 

The exposure length was adjusted by the brightness of 
each target, which was estimated during the short pointing 
exposures, depending also on the sky conditions. The aver¬ 
age photometric error in magnitudes is 0.04 mag. The data 
were processed using the specialized pipeline described in 


^ http://skinakas.physics.uoc.gr 


detail by King et al. (2014) along with the telescope control 
system. 

Since we have introduced a Galactic latitude cut select¬ 
ing objects with |6| > 10°, the average colour excess in the 
directions of our targets is relatively low, E{B — V) = 0.11"“ 
(Schlafly & Finkbeiner 2011), implying that the interstellar 
polarization is less than 1.0% on average (Serkowski et al. 
1975). The statistical uncertainty in the degree of polariza¬ 
tion is less than 1% in most cases, while the EVPA is typi¬ 
cally determined with a precision of 1° - 10° depending on 
the source brightness and fractional polarization. Detailed 
description of the instrument model and error analysis is 
given in King et al. (2014). 

In order to resolve the 180° ambiguity of the EVPA we 
followed a standard procedure (see e.g. Abdo et al. 2010a; 
Ikejiri et al. 2011; Kiehlmann et al. 2013), which is based 
on the assumption that temporal variations of the EVPA 
are smooth and gradual, hence adopting minimal changes of 
the EVPA between consecutive measurements. We define the 
EVPA variation as A0n = |6»n-|-l-6'n|-x/0-(6'n + l)2 +Cr(6'n)2, 
where 0n+i and 9^ are the n +1 and n-th points of the EVPA 
curve and g{6^+\) and cr(0n) are the corresponding errors of 
the position angles. If A^n > 90°, we shift the angle ^n+i 
by ± n X 180°, where the integer ± n is chosen in such a 
way that it minimizes A6n- If A^n < 90°, we leave ^n+i 
unchanged. 

Our first period of regular photometric and polarimetric 
monitoring of blazars started in 2013 July and lasted until 
the end of 2013 November. During the five-month period 
we obtained more than 1100 measurements of 101 objects 
from our sample almost uniformly spread over the observing 
season of each object. The median cadence and total season 
length for objects with At-median smaller than 20 d (includ¬ 
ing the June survey data, Pavlidou et al. 2014) is presented 
in Fig. 1, which is discussed in more detail in Sec. 3.3. 

2.3 Gamma-ray observations 

The gamma-ray data were obtained with the Large Area 
Telescope (LAT) onboard the Fermi gamma-ray space obser¬ 
vatory, which observes the entire sky every 3 h at energies of 
20 MeV - 300 GeV (Atwood et al. 2009). We analysed LAT 
data in the energy range 100 MeV < E < 100 GeV using the 
unbinned likelihood analysis of the standard Fermi analy¬ 
sis software package Science Tools v9r33p0 and the instru¬ 
ment response function P7REP_SOURCE_V15. Source 
class photons (evclass=2) were selected within a 15° re¬ 
gion of interest centred on a blazar. Cuts on the satel¬ 
lite zenith angle (< 100°) and rocking angle (< 52°) were 
used to exclude the Earth limb background. The diffuse 
emission from the Galaxy was modelled using the spatial 
model glU,em_vQ7>jrevl. The extragalactic diffuse and resid¬ 
ual instrumental backgrounds were included in the fit as an 
isotropic spectral template iso_source-vQ7>. The background 
models^ include all sources from the 2FGL catalogue within 
15° of the blazar. Photon fluxes of sources beyond 10° from 
the blazar and spectral shapes of all targets were fixed to 
their values reported in 2FGL. The source is considered to 

^ http://fermi.gsfc.nasa.gov/ssc/data/access/lat/ 
2yr_catalog/gll_psc_v07.xml 
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Table 2. Observational data for EVPA rotations detected by RoboPol in 2013. Columns (1),(2) - blazar identifiers; (3) - redshift; (4) 
- observational season length; (5) - average rotation rate; (6) - total amplitude of EVPA change; (7) - number of observations during 
rotation; (8) - time duration of the rotation; (9) - TeV emission flag according to TeVCat^"^ (“Y” means that the blazar has been detected 
in gamma rays with E > I TeV, “N” - otherwise); (10) - blazar subclass (LBL, IBL, HBL denote low, intermediate and high synchrotron 
peaked BL Lacertae objects, ESRQ - flat spectrum radio quasar). 


Blazar ID 

Survey 

name 

z 

Tobs 

(d) 

' At > 
(deg/d) 

^^max 

(deg) 

^ points 

Trot 

(d) 

TeV 

Class 

RBPLJ0136+4751 

OC 457 

0.8591 

59 

-6.6 

-225 

6 

34 

N 

FSRQi® 

RBPLJ0259+0747 

PKS0256-I-075 

0.893= 

72 

-4.8 

-180 

6 

38 

N 

FSRQi® 

RBPLJ0721-I-7120* 

S5 0716-1-71 

0.31® 

88 

-14.8 

-208 

11 

14 

Y 

LBLio 

RBPLJ0854-I-2006* 

OJ 287 

0.306'^ 

51 

-6.7 

-154 

10 

23 

N 

LBLi® 

RBPLJ1048-I-7143 

S5 1044-1-71 

1.15® 

142 

-9.0 

-188 

22 

21 

N 

- 

RBPLJ1555-I-1111 

PG 1553-1-113 

- 

129 

5.6 

128 

8 

23 

Y 

HBLio 

RBPLJ1558-I-5625 

TXS1557-I-565 

0.3® 

137 

7.2 

222 

9 

31 

N 

IBL?ii 

RBPLJ1806-I-6949 

3C 371 

0.05’^ 

143 

-16.5 

-347 

7 

21 

N 

LBLii 

RBPLJ1806-I-6949 

— ft — 

— ft — 

— ft — 

13.3 

238 

5 

18 

N 

— ft — 

RBPLJ1927-I-6117 

S4 1926-1-61 

- 

135 

-4.4 

-105 

6 

24 

N 

LBLi® 

RBPLJ2202-I-4216 

BL Lac 

0.069® 

137 

-51.0 

-253 

5 

5 

Y 

LBLi® 

RBPLJ2232-I-1143 

GTA 102 

I.O 37 I 

140 

-15.6 

-312 

8 

20 

N 

FSRQi® 

RBPLJ2232-I-1143 

— ft — 

— ft — 

— ft — 

-11.8 

-140 

6 

12 

N 

-ft- 

RBPLJ2243-I-2021 

RGB J2243-I-203 

- 

169 

-5.9 

-183 

5 

31 

N 

LBLi= 

RBPLJ2253-I-1608 

3C 454.3 

O. 859 I 

159 

-18.3 

-129 

4 

7 

N 

FSRQi® 

RBPLJ2311-I-3425 

B2 2308-1-34 

1.817® 

36 

3.3 

74 

20 

23 

N 

FSRQi® 


* Source belongs to sample (iii) 

^(Hewitt & Burbidge 1987);^(Murphy et al. 1993);®(Nilsson et al. 2008);^(Nilsson et al. 2010); 

®(Polatidis et al. 1995);®(Falco et al. 1998);^(de Grijp et al. 1992);®(Vermeulen et al. 1995); 

®(Wills & Wills 1976);^° (Donato et al. 2001);^^(Ghisellini et al. 2011);^^(Nieppola et al. 2006); 

^®(Fan et al. 2012);^^http://tevcat.uchicago. edu 

be detected if the test statistic, TS, provided by the analysis 
exceeds 10, which corresponds to approximately a Scr detec¬ 
tion level (Nolan et al. 2012). The systematic uncertainties 
in the effective LAT area do not exceed 10 per cent in the 
energy range we use (Ackermann et al. 2012). This makes 
them insignificant with respect to the statistical errors, that 
dominate over the short time-scales analysed in this paper. 

Moreover our analysis is based on the relative flux varia¬ 
tions. Therefore the systematic uncertainties were not taken 
into account. 

Different time bins tint, from 1 week to 25 d were used, 
depending on the flux density of the object. In order to make 
the analysis more robust we increased sampling of the pho¬ 
ton flux curves shifting centres of the time bins by tint/4 
interval from each other. This prevents losses of possible 
short-term events in the light curves and reduces the depen¬ 
dence of results on the particular position of the time bins. 

The oversampling introduces an autocorrelation in the pho¬ 
ton flux curves, which is however inessential for the analysis 
used in this work. 


3 RESULTS 

3.1 Detected rotations of EVPA 

The optical emission polarization plane of blazars is often 
variable even within the course of a single night. There is no 
objective physical definition of an EVPA rotation. Strictly 
speaking, any change of the EVPA between two measure¬ 
ments constitutes a rotation. However typically only high- 
amplltude (> 90°), smooth and well tracked variations of 
the EVPA are considered as rotations in the literature. 


We accept a swing between two consecutive EVPA 
measurements A9 = |0i+i — 6i\ as significant if A6 > 
^CT(6»i+l)2 -b a{9i)'^. We define as an EVPA rotation any 
continuous change of the EVPA curve with a total ampli¬ 
tude A^max > 90°, which is comprised by at least four mea¬ 
surements with significant swings between them. Start and 
end points of a rotation event are defined by a change of the 
EVPA curve slope A9i/Ati by a factor of 5 or a change of its 
sign. This definition is rather conservative, and is in general 
consistent with rotations reported in the literature. 

Using this definition, we identified 14 rotations of the 
EVPA in 12 blazars from the main sample during the sea¬ 
son of 2013 (see Table 2). This number is comparable to the 
number of previously known events of this type. Two more 
blazars with detected rotations, namely RBPLJ0721-I-7120 
and RBPLJ0854-I-2006, belong to the additional sample of 
hand-picked sources. These blazars/events were not included 
in the statistical or frequency analysis of the following sec¬ 
tions in this paper. The full season EVPA curves along 
with the evolution of the polarization degree and the R- 
band flux density, for all 14 blazars with detected rotations, 
are shown in Fig. 2 and listed in Table 2. The EVPA rota¬ 
tions are marked by filled black points. Clearly the events 
we have considered as rotations based on our criteria are the 
largest A^max rotation events that appear in these data sets. 
They are all characterized by smooth variations with a well- 
defined trend. Two events plotted in Fig. 2 do not follow the 
definition strictly. These are the rotation events detected in 
the data sets of RBPLJ1048+7143 and RBPLJ2311+3425. 
In both cases the rotations were interrupted by short, low 
amplitude, albeit significant swings in the opposite direction 
with respect to the overall rotation. Since both events are 
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Figure 2. Evolution of polarization degree, polarization position angle and /?-band magnitude for blazars with a detected rotation in 
the first RoboPol season. Periods of rotations are marked by filled black points. 


MNRAS 000, 1-16 (2015) 


R-band magnitude R-band magnitude R-band magnitude 
































6 D. Blinov et al 


RBPLJ1558 + 5625 


150p 

100- 


; -50- 
-100- 
-150- 

' 0,50 - 

: ^ 

• 0.45- 




f 


++ 


500 520 540 560 580 

JD (-2456000) 




RBPLJ1806+6949 

tB 

-e- 

-e- 

4> 

+ 


-100 

-150 

-200 

-250 

-300 

-350 


16.9 
17.oj 
17.1 £ 
17.2.0 
17.3 



■g -50 
2-100 
-150 
-200 


16.3 

16.4-d 

16.5 I 

16.6 E 
I6.7I 
16. 
16.9 


- 14.5 

t 15.0 

15.5 

16.0 


RBPLJ2202+4216 



o 



© 


© 

© 

© 

o \ 

© 

® © 

o 

© © 

a © ® © 

© © 

© 

© ' ® 

«> 0, 
® a.* ® 

© *• © 





• 


G 

O 

© 

© 

a® 

• 

Q3> ^ 

^ ©• 




• o 

Q 

/% 



- 13.5 

- 14.0 


S|- 1 - 

»- X 

0-^- 


t 

• 

• 

-t- 

% 

e , 

1 A >t> 

T ©5) © O @ 

1- 4 + '^’$0 

^ A 

s 1 

3 f ' 

7? $ 

6- *<1. 

5- ® O 


500 550 

JD (-2456000) 


16.4 

16.6 

16.8 

17.0 

17.2 


Figure 2 - continued (Continued) Evolution of polarization degree, polarization position angle and i?-band magnitude for blazars with 
a detected rotation in the first RoboPol season. Periods of rotations are marked by filled black points. 
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Figure 2 - continued (Continued) Evolution of polarization degree, polarization position angle and ij-band magnitude for blazars with 
a detected rotation in the first RoboPol season. Periods of rotations are marked by filled black points. 


well sampled these small deviations do not introduce any 
significant difference in the overall EVPA trend. Hence both 
events can be considered as single, large A^max rotations. In 
addition the RBPLJ2311+3425 event has an amplitude of 
~ 74°, which is less than the lower limit we accepted. How¬ 
ever the start and end points of the rotation are not defined 
due to a sparse sampling. It is likely that this well defined 
EVPA change would meet the 90° limit if we had a longer 
data set for this object. It is for this reason that we include 
this event in our sample of rotations. Both events have not 
been used in any of our statistical analyses involving com¬ 
parison between simulated and observed rotations. 

Some of the EVPA rotation events are coincident with 
an increase in the total flux, as it follows from a visual in¬ 
spection of Fig. 2. A quantitative comparison between the 
optical flux and the polarization variations will be presented 
in a forthcoming paper. 

3.2 General properties of EVPA rotations and 
rotators 

We estimated the maximal amplitude A^max and the dura¬ 
tion of the rotations Trot, using the first and last points of 
each event. Due to a moderate sampling and 180° EVPA am¬ 
biguity, the rotation start and/or end points cannot be pin¬ 
pointed accurately in five cases (namely RBPLJ0136-I-4761, 
RBPLJ0259+0747, RBPLJ1048+7143, RBPLJ1806+6949 
and RBPLJ2311-I-3425). This ambiguity affects the esti¬ 
mated A^max and Trot of the event, which should really be 
considered as lower limits in this case. We also estimated 
the average rotation rate as {A6/At) — A^max/Trot. These 
parameters as well as the blazar class and the TeV emission 
flag are listed in Table 2. 

We also collected data from the literature on previously 
known rotations of EVPA in blazars which show this be¬ 
haviour (“rotators” hereafter). Rates and A0max of these ro¬ 
tations were estimated from plots in the respective papers. 
These parameters as well as the blazar class and the TeV 
emission flag are listed in Table 3. 



Figure 3. Distributions of amplitudes and rates of EVPA ro¬ 
tations detected in RoboPol’s first season and reported in the 
literature. 


The distribution of A^max and rates of EVPA rotations 
from historical and RoboPol data are shown in Fig. 3. The 
number of detected rotations clearly decreases with grow¬ 
ing A^max. At the same time slow rotations dominate in the 
sample. This is presumably caused by a selection effect, be¬ 
cause fast rotations require better sampling of observations. 

Summarizing data on all known EVPA rotations in 
blazars to date we can list the following properties: 

(i) all known blazars with detected EVPA rotations are 
in the 2FGL catalogue (i.e. they are “gamma-ray-loud” 
sources); 

(ii) there are blazars known as TeV emitters as well as 
non-TeV sources among rotators; 

(iii) all subclasses of blazars show rotations of the EVPA, 
regardless of the position of the synchrotron peak maximum 
or the BL Lac/FSRQ dichotomy; 

(iv) there are eight blazars with more than one rotation 
detected. Comparison of these rotations shows that a single 
source can show rotations in both directions (five blazars 
known so far with this behaviour) and rotations observed 
in the same source can be of significantly different rates 
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Table 3. Data on known rotation of optical EVPA in blazars. Columns (1),(2) - blazar identifiers; (3) - 
average rotation rate; (4) - total amplitude of EVPA change; (5) - TeV emission flag according to TeVCat® 
(“Y” means that the blazar has been detected in gamma rays with E > 1 TeV, “N” - otherwise); (6) - blazar 
subclass (LBL, IBL, HBL denote low, Intermediate and high synchrotron peaked BL Lacertae objects, FSRQ 
- flat spectrum radio quasar); (7) - reference. 


Blazar ID 

Survey 

name 

\ At 1 

(deg d-l) 

^^max 

(deg) 

TeV 

Class 

Reference 

RBPLJ0423-0120 

PKS 0420-014 

-11.1 

-110 

N 

FSRQ4 

(D’Arcangelo et al. 2007) 

RBPLJ0721+7120 

S5 0716-1-71 

-1-130 

+ 180 

Y 

LBLi 

(Larionov et al. 2013b) 

RBPLJ0854+2006 

OJ 287 

-17 

-180 

N 

LBLi 

(Kikuchi et al. 1988) 

RBPLJ0958+6533 

S4 0954-1-65 

-1-18.2 

+240 

N 

LBL2 

(Larionov et al. 2011) 

RBPLJ1221+2813 

W Comae 

> -1-3.0 

+110 

Y 

IBL3 

(Benitez et al. 2013) 

RBPLJ1256-0547 

3C 279 

-9 

-180 

Y 

FSRQ-* 

(Abdo et al. 2010a) 

RBPLJ1256-0547 

3C 279 

-k4.3 

+290 

—// — 

— II — 

(Larionov et al. 2008) 

RBPLJ1256-0547 

3C 279 

-1-4.7 

+140 

—// — 

—// — 

(Aleksic et al. 2014a) 

RBPLJ1512-0905 

PKS 1510-089 

+15.6 

+720 

Y 

FSRQ'‘ 

(Marscher et al. 2010) 

RBPLJ1512-0905 

PKS 1510-089 

+12 

+400 

-II — 

— II — 

(Aleksic et al. 2014b) 

RBPLJ1512-0905 

PKS 1510-089 

-50 

-250 

-II — 

— II — 

(Aleksic et al. 2014b) 

RBPLJ1512-0905 

PKS 1510-089 

+11.7 

+500 

—n— 

—// — 

(Sasada et al. 2011)® 

RBPLJ2202+4216 

BL Lac 

+46 

+220 

Y 

IBLi 

(Marscher et al. 2008) 

RBPLJ2202+4216 

BL Lac 

+21 

+210 

—// — 

— II — 

(Sillanpaa et al. 1993) 

RBPLJ2232+1143 

CTA 102 

-60 

-180 

N 

FSRQ-* 

(Larionov et al. 2013a) 

RBPLJ2253+1608 

3C 454.3 

+16.3 

+130 

N 

FSRQ4 

(Sasada et al. 2010) 

RBPLJ2253+1608 

3C 454.3 

+9.3 

+400 

— // — 

— // — 

(Sasada et al. 2012) 


^(Donato et al. 2001);^(Ghiselllni et al. 2011);®(Tagliaferri et al. 2000);‘^(Fan et al. 2012) 
®http;//tevcat.uchicago. edu;®same as in Marscher et al. (2010). 


(in seven blazars rates differ by a factor larger than two 
in speed). 

3.3 Observed frequency of EVPA rotations 

The efficiency of an EVPA rotation detection depends on 
the intrinsic rate of the rotation as well as the frequency 
and uniformity of the observing cadence. The ambiguity of 
the polarization position angle introduces an upper limit on 
the rotation rate that can be unequivocally detected with a 
given typical cadence of observations. Clearly, for a typical 
time interval between observations (At), no EVPA rotation 
with a rate higher than 90°/(At) can be observed. 

For each blazar in our sample we found the median time 
difference between successive observations At-median and 
the total observing season length (defined as the time differ¬ 
ence between the first and the last observations) Tobs- These 
quantities (for blazars observed with At-median < 20 days) 
are shown in Fig. 1. In the same figure, we also plot three 
lines which indicate the necessary At-median and Tobs for 
detection of EVPA rotations at rates < 10 (solid line), < 15 
(dashed line) and < 20 (dotted line) degrees per day. 

The leftmost vertical part of each line represents the 
shortest Tobs needed to detect a rotation of A^max = 90° 
at a given rotation rate. The inclined portion of each line 
is determined by our requirement on a rotation event to be 
comprised by a minimum of four points. Given this require¬ 
ment, as At-median increases, so does Tobs- An EVPA data 
set with At-median and Tobs on that line can allow detec¬ 
tion of EVPA rotations with 90° < A0max < 270°. The 
horizontal part indicates the maximum At-median allowed 
the detection of a rotation event under the requirement of 
/S.9 < 90° in EVPA between two consecutive points. 

We can now estimate the frequency with which EVPA 


rotations appear in blazars as follows. Out of the 14 detected 
rotations in blazars of the main sample, 8 have rates less 
than 10 deg d“^. There are also 41 main sample (“gamma- 
ray bright”) blazars that were observed with At-median and 
Tobs (see Table 4) within the region defined by the solid 
line in Fig. 1. The total observing length for these blazars is 
6432 d. Thereby we estimate the frequency of “slow” rota¬ 
tions (rate < 10 deg d~^) in the main sample sources as one 
rotation in ~ 800 days (6432 d / 8 rotations). Following the 
same reasoning we estimate average frequencies of rotations 
for blazars in the main sample with rates < 15 deg d~^ and 
< 20 deg d“^ as one rotation in ~ 490 d (4912/10) and 
~ 180 d (2363/13), respectively. 

3.4 EVPA variability in blazars of different 
samples 

In order to address the question whether “the EVPA vari¬ 
ability is different in objects where rotations were detected 
compared to the rest of the main sample and to the control 
sample” we collate all EVPA “swing” events and measure 
their A^max and rates. We define an EVPA “swing” as any 
continuous change of the EVPA curve, without a lower limit 
in its A^max or in the number of measurements. As before, 
start and end points of a swing event are defined by a change 
of the EVPA curve slope by a factor of 5 or a change of its 
sign. 

We identified all such events for all blazars of the main 
and control samples within the 10 deg d~^ “detection box” in 
Fig. 1, and measured their amplitude, A^max, and mean ro¬ 
tation rate. The cumulative distribution function (hereafter 
CDF; e.g. Wall & Jenkins 2012) of the EVPA swings A^max 
and rotation rates for blazars in the main sample which 
showed rotations (“rotators”), blazars in the main sample, 
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Table 4. Sources of the main and control samples within the (■^) < 10 deg d ^ “detection box”. 


Blazar ID 
(RBPL. ..) 

Survey 

name 

Fobs 

(d) 

(At) 

(d) 

Blazar ID 
(RBPL.. .) 

Survey 

name 

Fobs 

(d) 

(At) 

(d) 


Main sample 



J1838+4802 

GB6 J1838+4802 

121 

7.0 

J0045+2127 

GB6 J0045-I-2127 

33 

4.0 

J1841+3218 

RXJ1841.7+3218 

152 

6.0 

J0114+1325 

GB6 J0114-I-1325 

38 

6.0 

J1903+5540 

TXS 1902+556 

135 

5.0 

J0211+1051 

MGl J021114-I-1051 

85 

4.0 

J1959+6508 

lES 1959+650 

143 

5.0 

J0217+0837 

ZS0214+083 

85 

6.0 

J2005+7752 

S5 2007+77 

140 

7.0 

J0423-0120 

PKS 0420-01 

12 

4.0 

J2015-0137 

PKS2012-017 

155 

6.5 

J0841+7053 

4C 71.07 

71 

6.0 

J2016-0903 

PMNJ2016-0903 

155 

7.0 

J1512-0905 

PKS 1510-08 

88 

2.0 

J2022+7611 

S5 2023+760 

158 

7.0 

J1542+6129 

GB6 J1542-I-6129 

87 

4.0 

J2030-0622 

TXS 2027-065 

143 

5.0 

J1553+1256 

PKS 1551+130 

132 

4.0 

J2039-1046 

TXS 2036-109 

144 

5.5 

J1604+5714 

GB6 J1604+5714 

135 

7.0 

J2131-0915 

RBS1752 

127 

5.0 

J1607+1551 

4C 15.54 

136 

8.0 

J2143+1743 

OX 169 

119 

5.0 

J1635+3808 

4C 38.41 

121 

2.0 

J2148+0657 

4C 6.69 

152 

4.5 

J1642+3948 

3G 345 

148 

6.0 

J2149+0322 

PKSB 2147+031 

169 

6.5 

J1653+3945 

Mkn 501 

153 

4.0 

J2150-1410 

TXS 2147-144 

130 

8.0 

J1725+1152 

IH 1720+117 

120 

3.0 

J2225-0457 

3C 446 

144 

6.0 

J1748+7005 

S4 1749+70 

87 

3.0 

J2251+4030 

MG4 J225201+4030 

177 

6.5 

J1751+0939 

OT 081 

154 

4.0 

J2334+0736 

TXS 2331+073 

138 

8.0 

J1754+3212 

RXJ1754.1+3212 

134 

5.0 

J2340+8015 

BZBJ2340+8015 

113 

5.5 

J1800+7828 

S5 1803+784 

133 

5.0 


Control sample 



J1809+2041 

RXJ1809.3+2041 

152 

4.5 

J1551+5806 

SBS1550+582 

118 

5.0 

J1813+3144 

B2 1811+31 

150 

6.0 

J1638+5720 

S4 1637+57 

138 

4.5 

J1836+3136 

RXJ1836.2+3136 

151 

6.0 

J2042+7508 

4G +74.26 

99 

5.0 


which did not show rotations (“non-rotators”), as well as for 
blazars in the control sample, are shown in Fig. 4. 

We performed a two sample Kolmogorov-Smirnov (K- 
S) test (e.g. Wall & Jenkins 2012) pairwise for three sam¬ 
ples of collected swing amplitudes and rates with the null- 
hypothesis that these samples are drawn from the same dis¬ 
tribution. The null-hypothesis is rejected for rotators and 
non-rotators with the p-value = 1.2 x 10“®, and for rota¬ 
tors and the control sample (p-value = 5 x 10“®). At the 
same time the distribution of swing amplitudes in the non¬ 
rotators and control sample sources is indistinguishable ac¬ 
cording to the test (p-value = 0.35). The maximum differ¬ 
ence between the CDFs of non-rotators and rotators is 0.29. 
It is reached at A^max ~ 25°. Even if we exclude the 14 
rotations (i.e. the largest A0max swings) of the main sample 
blazars, rotators still remain different from the non-rotators 
(p-value = 2 X 10“®). 

A similar analysis (as the one for A^max) for the distri¬ 
butions of EVPA swing rates leads to the same conclusion. 
The null-hypothesis is rejected for the rotators and the non¬ 
rotators (p-value = 1.4 x 10“®) and rotators vs. the control 
sample (p-value = 5 x 10“®), while it can not be rejected for 
the non-rotators and control sample (p-value = 0.18). 

We therefore conclude that blazars with detected rota¬ 
tions show significantly larger A^max and faster EVPA vari¬ 
ations when compared to blazars with no detected rotations. 
This difference cannot be attributed to differences in the 
sampling properties of the data sets. Therefore, the lack of 
detection of EVPA rotations in the “non-rotators” member 
of the main sample, as well as the blazar in the control sam¬ 
ple, may have a physical origin. Most of the non-rotators 
in the main and control samples may never show an EVPA 
rotation. 


4 RANDOM WALKS AS THE ORIGIN OF 

EVPA ROTATIONS 

4.1 MC simulations of EVPA swings 

Potentially EVPA swings can be explained by a stochastic 
process, which is physically justified by a presence of many 
independent cells in the emission region (e.g. Jones et al. 
1985; D’Arcangelo et al. 2007). According to this interpre¬ 
tation, the magnetic field is turbulent and apparent rotations 
result from a random walk of the full polarization vector di¬ 
rection as new cells with random magnetic field orientations 
appear in the emission region (Marscher 2014). In order to 
estimate the probabilities that the EVPA rotations we ob¬ 
served with RoboPol are produced by this kind of multicell 
random walk process we performed MC simulations of the 
stochastic variability of the polarization vector on the QU 
plane following Kiehlmann et al. (2013). 

For each blazar where an EVPA rotation event was 
observed, we created 10^ artificial light curves, each one 
with duration Fobs- The time steps Ati between consecutive 
points were drawn from a truncated power-law distribution, 
which approximates well the distribution of the time steps in 
all observed lightcurves. The parameters of this distribution 
(Atmin, Atmax and the power-law index) were determined by 
fitting it to the distribution of observed Ati for each object. 

The total flux density h emitted at each time step Ati, 
was drawn from a log-normal distribution. Such a distribu¬ 
tion approximates reasonably well the distribution of the ob¬ 
served flux densities for all blazars. The mean and variance 
of the log-normal distribution was set equal to the sample 
mean and variance of the distribution of the flux density of 
each blazar. 

The maximum possible fractional polarization produced 
by a uniform magnetic field is Pmax = (a -I- l)/(a -b 5/3) « 
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^^max [deg] {A8/At) [deg/day] 


Figure 4. CDFs of A6max and average rates for the main sample 
rotators (95 EVPA swings), the main sample non-rotators (298 
swings) and control sample sources (11 swings). See Sec. 3.4 for 
details. 


0.78 (Pacholczyk 1970). In the case of unresolved emis¬ 
sion region comprising N independent cells with a uni¬ 
form magnetic field, but randomly oriented among them, 
the average fractional polarization is given by the equation 
(Hughes & Miller 1991): 

(J’obs) « (1) 

We used this equation and the observed average polarization 
fraction, (Pobs), to estimate the number of cells, N, for each 
blazar. Each fc-th cell at i-th time step was assigned a flux 
density 7i,k (which was set equal to li/N for all cells at each 
time step) and a set of fractional Stokes parameters qi,k and 
iii^k- They were found as 


<li,k = <7i,k 




M.k = Ui,k 


( 2 ) 




where g^k and uf k are two numbers drawn from the standard 
normal distribution. Thereby the emission of each cell has 
polarization fraction Pmax- The sums Qi = p 9i,k and 
Ui = li determine the total Stokes parameters of 

the emitting region at each time step. 

At each time step the Stokes parameters of Avar(Ati) 
cells, selected randomly, were replaced by new values. The 
number of cells for replacement was estimated (from the 
average variance of the polarization degree) as follows: 


A'var(Ati) 


Ati CT (Pobs) 
(At) (Pobs) ’ 


(3) 


where (T(Pobs) is the observed standard deviation of the de¬ 
gree of polarization for each blazar, and (At) is the average 
time difference between observations. 

It was confirmed that the simulated and observational 
data in corresponding blazars have similar statistical prop¬ 
erties. Namely, the standard deviation and average of the 
polarization fraction are consistent with a{Pobs) and (Pobs). 


4-1.1 Individual rotations 

Using the algorithm described in Sec. 3, i.e. the same algo¬ 
rithm we used to identify rotations in real data, we iden¬ 
tified all rotations in the simulated data and found the 
number Wot of “successful” data sets, where at least one 


Table 5. Random walk modelling results for EVPA rotations 
detected by RoboPol in 2013. (1) - blazar identifier; (2) - occur¬ 
rence of rotations with AS^jax simul AS^jax obs estimated from 
the simulations; (3) - probability that a rotation produced by the 
random walk will be observed in Tobs ■ 


Blazar ID 

Tocc 

(d) 

P(RW) 

RBPLJ0136-I-4751 

505 

0.11 

RBPLJ0259-I-0747 

151 

0.48 

RBPLJ0721-I-7120 

325 

0.28 

RBPLJ0854-I-2006 

142 

0.36 

RBPLJ1048-I-7143 

180 

0.79 

RBPLJ1555-I-1111 

128 

1.00 

RBPLJ1558-I-5625 

266 

0.51 

RBPLJ1806-I-6949 

965 

0.15 

RBPLJ1806-I-6949 

259 

0.55 

RBPLJ1927-I-6117 

137 

0.98 

RBPLJ2202-I-4216 

633 

0.21 

RBPLJ2232-I-1143 

1557 

0.09 

RBPLJ2232-I-1143 

178 

0.87 

RBPLJ2243-I-2021 

183 

0.92 

RBPLJ2253-I-1608 

184 

0.86 

RBPLJ2311-I-3425 

61 

0.74 


rotation with A^max larger or equal to ASmax,obs was de¬ 
tected. We then estimated two ratios: P(RW) = Wot/10"* 
and Toco = 10"* • Pobs/Wot- The first ratio determines the 
probability to observe an EVPA rotation due to a random 
walk for each one of the observed EVPA curves for the given 
Tobs- The second ratio determines the average time interval 
between random walk rotations (i.e. the average occurrence 
rate for each blazar). The probabilities P(RW) and Tocc are 
listed in Table 5. The probabilities are larger than 10% in 
all but one object, and in some cases, they approach unity. 
This result indicates that the rotations we observed in some 
objects could be the result of a random walk process. 


4.1.2 Rotations as a population 

In this section we test the hypothesis that all the rotations 
observed by RoboPol in blazars of the main sample are pro¬ 
duced by the stochastic process. According to the analysis 
in Sec. 3.4 blazars exhibiting rotations have different prop¬ 
erties when compared to non-rotators. Therefore the sample 
of rotators must be considered separately. 

We performed the following simulation. At each itera¬ 
tion, an artificial EVPA curve was generated individually for 
each rotator from the main sample as explained in Sec. 4.1. 
In each of the simuated EVPA curves we identified the 
largest rotation and constructed the CDF of ASmax,simul 
among the blazars. An iteration was considered to be “suc¬ 
cessful” only in the case when the CDF of A^max,simul was 
lower or equal to the CDF of A0max,obs, i.e. the simulated set 
of EVPA curves had higher or equal fraction of rotations of 
a given length compared to the observed set. In the cases of 
RBPLJ1806-I-6949 and RBPLJ2232-I-1143 where double ro¬ 
tations were observed, we simulated only the largest A0max 
rotations. 

The CDF of A^max.obs along with a subset of 100 simu¬ 
lated CDFs is shown in Fig. 5. It was found that only 1.5% in 
10"^ trials were “successful”. Therefore, the probability that 
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Figure 5. CDFs of A^max in observed and a subset of 100 sim¬ 
ulated rotations. 

10 largest rotations in blazars of the main sample observed 
in our monitoring campaign all together were produced by a 
random walk is ~ 1.5%. If we repeat this simulation for the 
whole set of 16 EVPA rotations this probability is reduced 
to 0.5%. 

We conclude that, although some of the rotation events 
that we have detected may have been caused by a random 
walk process (as it is modelled in this paper), this hypothesis 
is not a likely explanation of the total number of detected 
EVPA rotations in our data set. 


5 OPTICAL EVPA ROTATIONS AND 

GAMMA-RAY ACTIVITY 

5.1 Average gamma-ray flux during EVPA 
rotations 

It has been suggested (e.g. Abdo et al. 2010a; 
Marscher et al. 2010) that rotations of EVPA in opti¬ 
cal emission of blazars are physically related to gamma-ray 
flares. 

In order to quantify a possible connection between 
EVPA rotations and gamma-ray activity, we first compared 
the average gamma-ray photon fluxes for each blazar dur¬ 
ing rotation events with the rest of the RoboPol season 
where no rotation was detected. Fig. 6 shows the gamma-ray 
lightcurves of blazars with detected rotations of EVPA. The 
green (light) area indicates the first year RoboPol observa¬ 
tional season for each object and the pink (dark) area indi¬ 
cates the period of the detected rotation. The average pho¬ 
ton fluxes (listed in Table 6) were calculated using the time 
intervals corresponding to the rotating and non-rotating pe¬ 
riods as single time bins (or averaging fluxes for two/three 
non-rotating time bins in cases, where they are split by the 
rotations). The gamma-ray photon flux during a rotation 
was higher than the flux during the rest of the season at l-cr 
level only in four cases. The average difference between the 
photon flux during rotations and along the rest of the season 
is —0.3 ± 3.4 X 10“®phcm“^ s“^. Thus, we do not observe 
any significant systematic change of the average gamma-ray 
photon flux simultaneous with the EVPA rotations. 

However, a comparison of the mean flux levels during 
the rotation and over a relatively long period may not be the 
best way to search for a correlation between the gamma-ray 


Table 6. Gamma-ray photon flux level during rotations and 
throughout the rest of the RoboPol season. 


Photon flux {E > 100 MeV) 
Blazar ID (10“^phcm“^s“^) 

rotation no rotation 


RBPLJ0136-I-4751 

0.40 

± 

0.14 

0.59 

± 

0.16 

RBPLJ0259+0747 

1.27 

± 

0.21 

< 

0 . 

71 

RBPLJ0721+7120 

0.95 

± 

0.18 

0.84 

± 

0.11 

RBPLJ0854+2006 

0.33 

± 

0.16 

0.91 

± 

0.18 

RBPLJ1048+7143 

3.39 

± 

0.32 

2.12 

± 

0.11 

RBPLJ1555+1111 

0.51 

± 

0.11 

0.54 

± 

0.05 

RBPLJ1558+5625 

< 

0 . 

34 

0.21 

± 

0.05 

RBPLJ1806+6949 

0.35 

± 

0.15 

0.40 

± 

0.07 

RBPLJ1806+6949 

< 

0 . 

83 

0.40 

± 

0.07 

RBPLJ1927+6117 

0.29 

± 

0.13 

0.09 

± 

0.05 

RBPLJ2202+4216 

4.67 

± 

0.93 

3.29 

± 

0.21 

RBPLJ2232+1143 

3.82 

± 

0.32 

3.34 

± 

0.25 

RBPLJ2232+1143 

4.55 

± 

0.70 

3.34 

± 

0.25 

RBPLJ2243+2021 

0.11 

± 

0.06 

0.17 

± 

0.04 

RBPLJ2253+1608 

6.98 

± 

0.68 

8.82 

± 

0.19 

RBPLJ2311+3425 

2.09 

± 

0.27 

2.13 

± 

0.35 


activity and EVPA rotations. For instance in the cases of 
RBPLJ0721-b7120 and first rotation of RBPLJ2232+1143, 
rotations are clearly coincident with prominent flares, al¬ 
though the average gamma-ray photon fluxes are indistin¬ 
guishable, since the season comprises a number of flaring 
events with similar amplitude. Moreover, rotations of EVPA 
can either precede or follow gamma-ray flares according to 
various theoretical scenarios. It is therefore important to 
search for a possible correlation between EVPA rotations 
and gamma-ray flares. 

5.2 Time lags between flares and EVPA rotations 

In order to investigate this relation, we first identified all 
flares that happened in the gamma-ray light curves during 
the RoboPol observing season. 

We adopted a formal definition of a gamma-ray flare 
similar to the one proposed by Nalewajko (2013): “a flare is 
a contiguous period of time, associated with a given photon 
flux peak, during which the photon flux exceeds half of the 
peak value, and this lower limit is attained exactly twice - at 
the start and at the end of the flare”. However the definition 
was slightly changed because Nalewajko (2013) analysed a 
sample of the brightest flares ever detected by Fermi LAT, 
while we are interested in even smaller amplitude events. We 
found that a peak photon flux excess factor equal to 2/3, 
instead of the original 1/2 proposed by Nalewajko (2013), 
gives a better agreement with a visual identification of flares 
in the photon flux curves. Intervals of the photon flux curves 
identified as flares are marked by red (light) points in Fig. 6. 

We searched for the closest gamma-ray flare to the ro¬ 
tation event of each rotator, and we fitted it using a profile 
with an exponential rise and decay. This kind of profile is 
commonly used for fitting an individual blazar flare pulse 
in optical, gamma and radio bands (e.g. Abdo et al. 2010b; 
Chatterjee et al. 2012): 

/ tp—t t—tp \ ^ 

F{t) ^Fp + Fpl e^ -b e Ya- j , (4) 
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Figure 6. Gamma-ray light curves of objects with detected rotations of EVPA. The RoboPol observational season is marked by the 
green (light) area. The pink (dark) area shows duration of the rotation. Green ticks mark moments of our optical EVPA measurements. 
All curves are centred to the mean day of the RoboPol observing season. Detected flares are marked by red points, while the blue curve 
is the analytical function fit of the flares closest to observed rotations (see text for details). Vertical dashed lines indicate intervals of the 
light curves used in the fitting procedure. 
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Figure 7. Distributions of observed Tots s-nd simulated Tginmi 
time lags between middle points of rotations and tp of gamma- 
ray flares. 


where Fc represents an assumed constant level underlying 
the flare, Fp measures the amplitude of the flare, tp de¬ 
scribes the time of the peak (it corresponds to the actual 
maximum only for symmetric flares), and T'* measure the 
rise and decay time, respectively. All the parameters were set 
to be free, while initial values used in the fitting procedure 
were estimated from the photon flux curves. Upper limits 
of the light curve were not used in the fitting procedure. 
In the case of RBPLJ0721+7120 three flares that occurred 
during the observing season were fitted simultaneously be¬ 
cause a single flare fit resulted in an unrealistic F^ value. In 
the cases of double rotations in RBPLJ1806-I-6949 and RB- 
PLJ2232-I-1143 the flares closest to the rotations were also 
fitted together to provide a consistent Fc value. In addition, 
in three cases, the closest flares happened just outside the 
RoboPol season interval. The best fitting curves are shown 
in Fig. 6. 

We estimated time lags, robs, between rotations and 
the closest gamma-ray flares as robs = — tp, where 

T''°* is the middle point of each EVPA rotation, defined 
as trot,start + (see Table 7). The time lags have a dis¬ 
tribution, shown by green (light) bars in Fig.7, with mean 
and standard deviation equal to 5.1 and 21.8 d, respectively. 
The distribution is indistinguishable from the normal distri¬ 
bution Af(0, 21.8) following the K-S test (p-value = 0.39). 
Thereby we do not find any preference for positive or neg¬ 
ative robs. A distribution of observed time lags is expected 
to be close to a normal distribution with the mean at zero 
if rotations of the EVPA are not connected to gamma-ray 
flares. Because, in this case, the overall distribution is pro¬ 
duced by a set of random values each having distributions of 
different widths and symmetric with respect to zero. For this 
reason the time lags distribution does not, on its own, sup¬ 
port a physical connection between gamma-ray flares and 
rotations. However theoretical models allow for either posi¬ 
tive or negative lags, depending on conditions and emission 
region properties, when a physical connection between rota¬ 
tions and gamma-ray flares does exist (Zhang et al. 2014). 
Therefore, a physical connection cannot be excluded based 
on the distribution of Tobs- 


Table 7. Gamma-ray flares fitting results. (1) - blazar identifier; 
(2) - time difference, Tobs, between tp of the closest gamma-ray 
flare and middle point of the rotation (positive means leading 
flare); (3) - gamma-ray flare amplitude measured relative to the 
average photon flux of the blazar from 2FGL. 


Blazar ID 


RBPLJ0136-I-4751 
RBPLJ0259-I-0747 
RBPLJ0721-I-7120 
RBPLJ0854-I-2006 
RBPLJ1048-I-7143 
RBPLJ1555-I-1111 
RBPLJ1558-I-5625 
RBPLJ1806-1-6949 
RBPLJ1806-I-6949 
RBPLJ1927-I-6117 
RBPLJ2202-I-4216 
RBPLJ2232-I-1143 
RBPLJ2232-I-1143 
RBPLJ2243-I-2021 
RBPLJ2253-I-1608 
RBPLJ2311-1-3425 


Tobs 7 -flare 
(d) rel. ampl. 


53.8 

0.6 ±0.08 

-2.4 

15.1 ± 2.9 

0.8 

1.0 ±0.5 

-5.3 

2.5±1.1 

2.5 

7.3 ±3.6 

42.9 

1.1 ±0.2 

-14.8 

1.9 ±0.9 

5.4 

0.7 ±0.6 

-27.8 

1.3 ±0.4 

7.5 

0.6 ±0.2 

3.1 

3.1 ±0.6 

2.6 

7.3 ±1.5 

-3.7 

12.1 ± 1.5 

29.0 

0.7 ±0.4 

-30.2 

1.7±0.2 

18.8 

16.6 ± 1.3 


5.3 Relation of gamma-ray flare amplitudes and 
time delays 

We normalized the amplitude, Fp, of the gamma-ray flare 
closest to the EVPA rotation event by the average photon 
flux of each blazar (as listed in 2FGL; Nolan et al. 2012). 
Corresponding values are listed in Table 7 and plotted as 
a function of Tobs in Fig. 8. The filled black squares show 
redshift-corrected time lags, i.e. Tcorr = Tobs/(l + z), while 
open circles show Tobs for blazars with unknown z. The “er¬ 
rors” on the time lags are defined as fitting errors of tp plus 
the time difference between the first/last point of a rota¬ 
tion event, and the previous/next closest point of the EVPA 
curve. Due to the lack of data in some cases these “uncertain¬ 
ties” are undefined, while in others, due to sparse sampling, 
they are almost certainly overestimated. 

A noticeable feature is that Tcorr is in the range (—6, +6) 
d for the most prominent gamma-ray flares. Basically, all 
five brightest flares have happened almost simultaneously 
with EVPA rotation events. The brightest flare which has 
the largest deviation from the zero-delay is the flare of RB- 
PLJ2311-I-3425 where the start point of the rotation is un¬ 
defined and therefore the time-delay has a large uncertainty. 

There are three more flares with similarly small time 
lags, and small relative amplitudes. Thus a small separation 
between a flare and a rotation is not a sufficient condition 
for extraordinary brightness of the high-energy flare. 

Separating the flares into two subsamples of high and 
low amplitude events (dashed line in Fig. 8) we examined the 
significance of the difference in time delays between them. 
The mean of the absolute Tobs values for the high and low 
amplitude subsamples is 5.2 and 20.1 d, respectively. Ac¬ 
cording to the Student’s t-test (e.g. Wall & Jenkins 2012), 
the difference between the two mean values is somewhat sig¬ 
nificant (p-value = 0.025). 
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Figure 8. Time lags, ro^si versus normalized gamma-ray flare 
amplitude, _Fp. Redshift corrected and non-corrected r^bs values 
are plotted with filled squares and open circles, respectively. 


5.4 Accidental proximity of rotations and 
gamma-ray flares 

5 . 4.1 Individual blazars 

In order to estimate the probability of the accidental ob¬ 
served proximity in time of rotations and gamma-ray flares, 
we performed MC simulations using the observed gamma- 
ray photon flux curves. This allows us to account for the 
real variability of blazars in the gamma-ray band. For each 
rotator we processed a long-term set of Fermi LAT data 
(54683 < MJD < 57065) with time bins equal to the ones 
used in Sec. 5.2. Then we identified and fitted all gamma-ray 
flares following the procedure described previously, using the 
same photon flux excess factor of 1.5. The number of flares 
identified in the photon flux curves of rotators is in the range 
of 12 - 76. After that we randomly assigned the middle point 
of a simulated rotation to a time on the photon flux curve 
and measured the time lag between the rotation and the 
closest gamma-ray flare, Tsimui. Repeating this simulation 
10"* times for each blazar, we determined the distributions 
of time delays Tsimui. Using these distributions, we estimated 
the probability of robs to be produced by chance P(robs), by 
calculating the fraction of simulations where Tsimui < Tobs- 
The probabilities range between 3 and 78 per cent (see Ta¬ 
ble 8). Pink (dark) boxes in Fig. 7 indicate the distribution of 
"rsimui, using the results from the simulation for all blazars. 
According to the K-S test the null hypothesis that Tsimui 
and Tobs are drawn from the same distribution can not be 
rejected (p-value = 0.38). Therefore, it is possible that the 
Tobs values we observed, may be accidental for each of the 
blazars in the sample. 

In Sect. 4.1 we determined the probability of the EVPA 
rotations to be observed in our observing window assuming 
that they are produced by a stochastic process. The simula¬ 
tions described above give us the probability of an acciden¬ 
tal simultaneity between these rotations and gamma-flares. 
Therefore the probability of superposition of both indepen¬ 
dent events: (a) random rotation and (b) random proximity 
to a gamma-ray flare, can be estimated as a product of the 


Table 8. Modelling results for the connection between EVPA 
rotations detected by RoboPol in 2013 and gamma-ray flares. (1) 
- blazar identifier; (2) probability of an accidental time lag; (3) - 
combined probability of a rotation being produced by the random 
walk and located as close to the corresponding gamma-ray flare 
as it was observed. 


Blazar ID 

P(Tobs) 

P(RW-bTobs) 

RBPLJ0136-I-4751 

0.75 

0.08 

RBPLJ0259-I-0747 

0.03 

0.02 

RBPLJ0721-I-7120 

0.04 

0.01 

RBPLJ0854-I-2006 

0.23 

0.08 

RBPLJ1048-1-7143 

0.14 

0.11 

RBPLJ1555-I-1111 

0.72 

0.72 

RBPLJ1558-I-5625 

0.20 

0.10 

RBPLJ 1806-1-6949 

0.10 

0.02 

RBPLJ 1806-1-6949 

0.49 

0.27 

RBPLJ1927-I-6117 

0.08 

0.08 

RBPLJ2202-I-4216 

0.21 

0.04 

RBPLJ2232-I-1143 

0.14 

0.01 

RBPLJ2232-I-1143 

0.19 

0.17 

RBPLJ2243-I-2021 

0.48 

0.44 

RBPLJ2253-I-1608 

0.78 

0.67 

RBPLJ2311-1-3425 

0.56 

0.41 


respective probabilities. These combined probabilities are 
less than 5 per cent for five events (see column 3 of Table 8) . 
This result indicates that, at least for some rotations, the 
random walk model and the absence of any physical connec¬ 
tion between the EVPA variability and high-energy activity 
is an unfavourable interpretation. 


5 . 4.2 Rotators as a population 

In order to assess the probability that the entire set of the 
time lags appeared in the main sample rotators in a random 
way, we run the following simulation. Repeating the proce¬ 
dure described in Sec. 5.2 we identified and fitted all flares in 
the gamma-ray photon flux curve (54683 < MJD < 57065) 
of each blazar from the main sample with a detected rota¬ 
tion. Then placing a simulated rotation at a random position 
on each of the gamma-ray curves, we defined the shortest 
time lag between the central point of the rotation and tp 
of the nearest flare. After this the CDF of absolute values 
of the simulated time lags was constructed for the set of 14 
events. 

Repeating the routine 10® times we found that only one 
out of every 5000 simulations produces a CDF which is in 
its entirety located closer to zero or coincides with the CDF 
of observed time lags (see Fig. 9). Thereby we estimate the 
probability that all 14 delays together were produced by 
chance as 2 x 10“^. When we repeat this procedure for all 
16 rotations together including two non-main sample events, 
the estimated probability decreases to 5 x 10“®. Therefore, 
it is very unlikely that none of the observed EVPA rotations 
is related physically to the flaring activity in gamma-rays. 
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Figure 9. CDFs of the time lags between the EVPA rotations 
middle points and tp of the closest gamma-ray flares for the main 
sample rotators. Black line — observed time lags, thin grey lines - 
10 "* simulated values for the whole sample of rotations (see text 
for details). 


6 CONCLUSIONS 

During the first season of operation of the RoboPol project, 
we detected 16 rotations of the polarization plane in op¬ 
tical emission of blazars. These detections double the ex¬ 
isting list of such events. All EVPA rotations are observed 
in blazars which are detected by Fermi, in agreement with 
previous experiments, which have detected similar events in 
the same class of objects. Our strategy of monitoring both 
gamma-ray-loud and quiet samples, suggests that the lack 
of EVPA rotations detection by RoboPol in gamma-ray- 
quiet objects cannot be due to a difference in the sampling 
pattern. Combining our results with those reported in the 
literature we found that rotations can be detected in both 
TeV and non-TeV emitters. Our results also indicate that all 
subclasses of blazars show rotations of the EVPA (regardless 
of the position of the synchrotron peak maximum or the BL 
Lac/FSRQ dichotomy). We expect that the results after the 
3-yr planned RoboPol monitoring campaign will allow an 
accurate determination of the rotations rate in the various 
blazar subclasses. 

Analysis of the first-year data shows that blazars with 
detected rotations have significantly faster and longer EVPA 
swings when compared to non-rotators. This suggests that 
rotations of EVPA may be specific for a particular activity 
state or for a subclass of blazars with peculiar properties. 

The fact that EVPA rotations have been detected only 
in gamma-ray-loud objects already suggests a physical rela¬ 
tion between gamma-ray and optical polarization variability 
in blazars. Nevertheless, we used extensive MC simulations 
to investigate whether the EVPA rotations we observed can 
be produced by a random walk process of the polarization 
vector. We found that a random walk process can result in 
EVPA rotations with A^max.simui as large as A^max.obs for 
the given At-median and Tobs of the individual RoboPol 
data sets. However, we also found that it is unlikely (proba¬ 
bility is < 1.5 X 10~^), that all the rotations that we observed 
in the first RoboPol season are due to a random walk pro¬ 
cess. 


The average gamma-ray photon fluxes do not show any 
significant systematic increase during the rotation events. 
We also found that, the time lags between rotations of the 
EVPA and nearest gamma-ray flares follow a Gaussian dis¬ 
tribution with a mean ~ zero. 

We performed a second set of MC simulations in order 
to assess the randomness of the observed time delays. Our 
results suggest that, on an individual basis, the time lags we 
observe do not necessarily suggest a physical link between 
EVPA rotations and gamma-ray flares. On the other hand, 
when we consider the rotators as a population, it is highly 
unlikely {p = 2x 10“"^) that the proximity of EVPA rotations 
to gamma-ray flares is accidental in all cases. Therefore at 
least some EVPA rotations must be physically connected to 
the high-energy activity. 

Our data suggest that, the highest amplitude gamma- 
ray flares may be physically connected with EVPA rotations, 
based on the fact that they are associated with smaller- 
than-average time lags. Perhaps there are two different types 
of gamma-ray flares, produced by different physical mecha¬ 
nisms. One of them may result in higher (than average) am¬ 
plitude flares and EVPA rotation events. The other one may 
produce the rest of the smaller amplitude flares, which are 
not related with the remaining rotations, probably produced 
by a random walk process. 

For the first time we studied a set of EVPA rotations 
discovered in a large, well-defined, regularly monitored sam¬ 
ple of blazars. The diversity of results found in individual 
and population analysis shows the importance of these kinds 
of studies. The RoboPol monitoring of blazars will continue 
for at least two more years. The question about the mech¬ 
anisms responsible for the EVPA rotations in blazars and 
their possible connection to the high-energy activity will be 
explored in more detail after accumulation of a larger data 
set by RoboPol. 

The statistical analysis of our data set required us to 
make subjective choices regarding the details of our event 
definitions and the test statistics we used. These, for exam¬ 
ple, include the definition of a rotation, the definition of a 
“gamma-ray activity” (increase in gamma-ray flux during a 
rotation versus proximity to a gamma-ray flare peak), the 
definition and fitting procedure of a gamma-ray flare, use of 
the ASmax and |robs| CDFs as test statistics and so forth. 
Making these choices introduces unavoidable and unquan- 
tifiable biases in our final results. However, our exploratory 
analysis of the first-year data presented here has allowed us 
to identify well-defined statistical questions, which we can 
address in a robust, a priori fashion using our second- and 
third-year data. 
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